Setup

Set up workspace, set options, and load required packages.

rm(list=ls(all=TRUE)) 
knitr::opts_chunk$set(root.dir = "~/R",warning=FALSE, message=FALSE)
library(gridExtra)
library(multcomp)
library(emmeans)
library(ggplot2)
library(dplyr)
library(effects)
library(plyr)
library(lattice)
library(glmmTMB)
library(effects)
library(lsmeans)
library(MuMIn)
library(MASS)
library(car)
library(FSA)
library(lmerTest)
library(lme4)
library(blmeco)
library(ggsci)
library(coxme)
library(survival)
library(cowplot)
library(RColorBrewer)
library(tidyverse)
library(partitionBEFsp)
library(factoextra)
library(ggfortify)
library(cowplot)
library(survminer)

1. GROW OUT SURVIVORSHIP

Load data

Survivorship is measured as binary on each individual spat in the experiment.

#load data
rearing<-read.csv("Data/GrowOut_Responses.csv", header=TRUE, sep=",", na.strings="NA")
rearing$Quadrant<-as.factor(rearing$Quadrant)
rearing$Diversity.Type<-relevel(rearing$Diversity.Type, ref="Individual")

#rearing$Total.Juveniles<-as.factor(rearing$Total.Juveniles)
#rearing$Fusion.Partners<-as.factor(rearing$Fusion.Partners)
#rearing$Richness<-as.factor(rearing$Richness)

Subset the final timepoint for additional analyses.

final<-rearing[which(rearing$Timepoint == "Final"),]
multiple<-rearing[which(rearing$Community == "Multiple"),]
multiple_final<-multiple[which(multiple$Timepoint == "Final"),]

Analysis

Analyze with a logistic binomial regression model.

Analyze at the final timepoint.

#with final dataset
survmod2<-glmer(Survivorship ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Genotype) + (1|Tank) + (1|Slide), data=final, family=binomial)

#final dataset
summary(effect(c("Richness", "Fusion.Partners"), xlevels=4, survmod2))
## 
##  Richness*Fusion.Partners effect
##         Fusion.Partners
## Richness         0         1         2         3
##        1 0.4749728 0.9356204 0.9957349 0.9997334
##        2 0.5151037 0.9175013 0.9914845 0.9991803
##        3 0.5550408 0.8948559 0.9830701 0.9974823
##        4 0.5942795 0.8668961 0.9666214 0.9922938
## 
##  Lower 95 Percent Confidence Limits
##         Fusion.Partners
## Richness         0         1         2         3
##        1 0.2857530 0.7943532 0.9407520 0.9836550
##        2 0.3608520 0.8075419 0.9468818 0.9860349
##        3 0.3649140 0.7750230 0.9272927 0.9767284
##        4 0.3218885 0.6699355 0.8147196 0.8865651
## 
##  Upper 95 Percent Confidence Limits
##         Fusion.Partners
## Richness         0         1         2         3
##        1 0.6716629 0.9820395 0.9997088 0.9999957
##        2 0.6665306 0.9671889 0.9986868 0.9999525
##        3 0.7303142 0.9545994 0.9962318 0.9997327
##        4 0.8188353 0.9543351 0.9947841 0.9995289
summary(effect("Fusion.Partners", xlevels=4, survmod2))
## 
##  Fusion.Partners effect
## Fusion.Partners
##         0         1         2         3 
## 0.5167374 0.9166729 0.9912419 0.9991419 
## 
##  Lower 95 Percent Confidence Limits
## Fusion.Partners
##         0         1         2         3 
## 0.3625326 0.8073334 0.9467831 0.9859895 
## 
##  Upper 95 Percent Confidence Limits
## Fusion.Partners
##         0         1         2         3 
## 0.6678192 0.9665337 0.9986131 0.9999481
summary(effect("Community", xlevels=4, survmod2))
## 
##  Community effect
## Community
## Individual   Multiple 
##  0.9580409  0.7382699 
## 
##  Lower 95 Percent Confidence Limits
## Community
## Individual   Multiple 
##  0.8545984  0.5848937 
## 
##  Upper 95 Percent Confidence Limits
## Community
## Individual   Multiple 
##  0.9888517  0.8495526
dispersion_glmer(survmod2) #no overdispersion
## [1] 0.9350818
plot(residuals(survmod2)) + abline(h=0, lty=2) #Dispersed randomly

## integer(0)

Create plot for effect of community type.

eff.surv1 <- predictorEffect("Community", survmod2)
surv.effects1<-plot(eff.surv1,
                   lines=list(multiline=TRUE, 
                              col=c("black"), 
                              lty=1), 
                   confint=list(style="bars", width=4, col="black"), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="right",
                   xlab=expression(bold("Number of Juveniles")), 
                   main=""); surv.effects1

                   #lattice=list(key.args=list(space="right",
                                              #border=FALSE, 
                                              #title=expression(bold("Genotypic \nRichness")),
                                              #cex=1, 
                                              #cex.title=1)));surv.effects1

Create plots for fusion partner effects.

survmod2a<-glmer(Survivorship ~ Richness + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Genotype) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Multiple"))

eff.surv2 <- predictorEffect("Fusion.Partners", survmod2a, xlevels=4, rug=TRUE)
surv.effects2<-plot(eff.surv2,
                   lines=list(multiline=TRUE, 
                              col=c("lightblue","mediumblue", "blue4", "darkblue"), 
                              lty=1), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="right",
                   xlab=expression(bold("Fusion Partners")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Genotypic \nRichness")),
                                              cex=1, 
                                              cex.title=1)));surv.effects2

Generate mean values of survivorship for 1) number of juveniles and 2) within “multiple” juvenile groups, mean values for number of fusion partners.

meansurv <- ddply(final, c("Community"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meansurv
##    Community   N      mean        sd         se max     lower    upper
## 1 Individual  58 0.7758621 0.4206554 0.05523477   1 0.3552066 1.196518
## 2   Multiple 163 0.6687117 0.4721270 0.03697984   1 0.1965847 1.140839
meansurv2 <- ddply(multiple_final, c("Fusion.Partners"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meansurv2
##   Fusion.Partners  N      mean        sd         se max      lower    upper
## 1               0 75 0.3466667 0.4791133 0.05532324   1 -0.1324466 0.825780
## 2               1 48 0.9583333 0.2019409 0.02914766   1  0.7563924 1.160274
## 3               2 24 0.9166667 0.2823299 0.05763034   1  0.6343368 1.198997
## 4               3 16 0.9375000 0.2500000 0.06250000   1  0.6875000 1.187500

2. GROW OUT FUSION

Analysis

Analyze with a logistic binomial regression model.

Final timepoint dataset.

Analyze fusion as a product of genotypic richness within “multiple” juvenile groups.

#with final dataset
fusemod1<-glmer(Fusion ~ Richness + (1|Parent.Site/Genotype) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Multiple"))
summary(fusemod1)
Anova(fusemod1, type="II") # no effects on rate of fusion by the end point

dispersion_glmer(fusemod1) #no overdispersion
plot(residuals(fusemod1)) + abline(h=0, lty=2)

There is an effect of genotypic diversity on fusion rates. Plot the relationship.

eff.fuse1 <- predictorEffect("Richness", fusemod1, xlevels=4, rug=TRUE)
fuse.effects1<-plot(eff.fuse1,
                   lines=list(multiline=TRUE, 
                              col=c("black"), 
                              lty=1),
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   #axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Successful Fusion)")), 
                   xlab=expression(bold("Richness")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Genotypic \nRichness")),
                                              cex=1, 
                                              cex.title=1)));fuse.effects1

Generate a summary of proportion fused at the end of the grow-out period.

meanfusion <- ddply(multiple_final, c("Richness"), summarise,
                   N    = length(Fusion[!is.na(Fusion)]),
                   mean = mean(Fusion, na.rm=TRUE),
                   sd   = sd(Fusion, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Fusion, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanfusion
##   Richness  N      mean        sd         se max       lower     upper
## 1        1 57 0.5263158 0.5037454 0.06672270   1  0.02257042 1.0300612
## 2        2 26 0.6153846 0.4961389 0.09730085   1  0.11924568 1.1115236
## 3        3 36 0.3888889 0.4944132 0.08240221   1 -0.10552434 0.8833021
## 4        4 44 0.6363636 0.4866071 0.07335878   1  0.14975654 1.1229707

3. GROW OUT GROWTH

Analysis

Growth is measured as polyp expansion in each spat by number of spat in each sibling type and number of fused individuals. Growth is measured in final dataset as all spat started with 1 polyp each.

Analyze polyp expansion with the final dataset and present means.

polypmod1<-glmer(Polyps ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Genotype) + (1|Tank) + (1|Slide), data=final, family = poisson)
summary(polypmod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Polyps ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners +  
##     (1 | Parent.Site/Genotype) + (1 | Tank) + (1 | Slide)
##    Data: final
## 
##      AIC      BIC   logLik deviance df.resid 
##    588.8    616.3   -285.4    570.8      148 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.11373 -0.54677 -0.02775  0.43720  1.40877 
## 
## Random effects:
##  Groups               Name        Variance  Std.Dev. 
##  Slide                (Intercept) 0.000e+00 0.000e+00
##  Genotype:Parent.Site (Intercept) 7.465e-02 2.732e-01
##  Parent.Site          (Intercept) 1.629e-10 1.276e-05
##  Tank                 (Intercept) 9.270e-03 9.628e-02
## Number of obs: 157, groups:  
## Slide, 92; Genotype:Parent.Site, 20; Parent.Site, 6; Tank, 5
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               1.22325    0.14229   8.597   <2e-16 ***
## Richness                  0.01432    0.07830   0.183    0.855    
## CommunityMultiple        -0.10036    0.17493  -0.574    0.566    
## Fusion.Partners          -0.03007    0.16617  -0.181    0.856    
## Richness:Fusion.Partners  0.01653    0.04814   0.343    0.731    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Rchnss CmmntM Fsn.Pr
## Richness    -0.545                     
## CmmntyMltpl  0.073 -0.687              
## Fusn.Prtnrs -0.337  0.651 -0.723       
## Rchnss:Fs.P  0.408 -0.783  0.669 -0.934
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(polypmod1, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Polyps
##                           Chisq Df Pr(>Chisq)
## Richness                 0.5259  1     0.4683
## Community                0.3291  1     0.5662
## Fusion.Partners          0.1522  1     0.6964
## Richness:Fusion.Partners 0.1179  1     0.7314
qqPlot(residuals(polypmod1))

## 813 817 
## 144 145
meangrowth <- ddply(final, c("Community"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Polyps, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meangrowth
##    Community   N     mean       sd        se max    lower    upper
## 1 Individual  45 3.533333 1.778661 0.2651472   8 1.754672 5.311995
## 2   Multiple 112 3.357143 1.558800 0.1472928   7 1.798343 4.915943

Plot the model outputs.

eff.polyps1 <- predictorEffect("Fusion.Partners", polypmod1, xlevels=4, rug=TRUE)
polyp.effects1<-plot(eff.polyps1,
                   lines=list(multiline=TRUE, 
                              col=c("lightblue","mediumblue", "blue4", "darkblue"), 
                              lty=1), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   #axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Total Polyp Expansion")), 
                   legend.position="right",
                   xlab=expression(bold("Fusion Partners")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Genotypic \nRichness")),
                                              cex=1, 
                                              cex.title=1)));polyp.effects1

Combine figures for Grow-Out Period.

grow_out_figs<-plot_grid(surv.effects1, surv.effects2, labels = c("A", "B"), ncol=2, nrow=1, rel_heights= c(1,1), rel_widths = c(0.8,1), label_size = 20, label_y=1, align="h");grow_out_figs

ggsave(filename="Figures/grow_out_figs.png", plot=grow_out_figs, dpi=500, width=10, height=6, units="in")

4. STRESS SURVIVORSHIP

Load data

Survivorship is measured as binary on each individual spat in the experiment.

#load data
stress<-read.csv("Data/Stress_Responses.csv", header=TRUE, sep=",", na.strings="NA")
stress$Diversity.Type<-relevel(stress$Diversity.Type, ref="Individual")
stress$Polyps<-as.numeric(stress$Polyps)

Subset the final timepoint, “multiple” juvenile datasets for additional analyses.

final_stress<-stress[which(stress$Timepoint == "End"),]
multiple_stress<-stress[which(stress$Community == "Multiple"),]
multiple_final_stress<-multiple_stress[which(multiple_stress$Timepoint == "End"),]
d26_stress<-stress[which(stress$Timepoint == "S2D11"),]

Analysis

Analyze with a logistic binomial regression model.

Check for effect of community type (“individual” or “multiple” juveniles) on survivorship and display means.

#with full time dataset
survmod3b<-glmer(Survivorship ~ Day * Treatment * Community + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide), data=stress, family=binomial, link=logit) 

summary(survmod3b)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Survivorship ~ Day * Treatment * Community + (1 | Parent.Site/Genotype) +  
##     (1 | Treatment:Tank) + (1 | Slide)
##    Data: stress
## 
##      AIC      BIC   logLik deviance df.resid 
##    446.3    506.2   -211.2    422.3     1076 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -129.711    0.001    0.026    0.116    6.235 
## 
## Random effects:
##  Groups               Name        Variance  Std.Dev.
##  Slide                (Intercept) 2.523e+00 1.588500
##  Genotype:Parent.Site (Intercept) 5.680e-01 0.753637
##  Treatment:Tank       (Intercept) 3.662e-06 0.001914
##  Parent.Site          (Intercept) 8.331e-04 0.028863
## Number of obs: 1088, groups:  
## Slide, 90; Genotype:Parent.Site, 20; Treatment:Tank, 9; Parent.Site, 6
## 
## Fixed effects:
##                                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          6.57183    1.66342   3.951 7.79e-05 ***
## Day                                 -0.08040    0.05789  -1.389 0.164885    
## TreatmentHigh                        0.97861    1.88333   0.520 0.603332    
## CommunityMultiple                    3.65815    2.44562   1.496 0.134708    
## Day:TreatmentHigh                   -0.28864    0.08140  -3.546 0.000391 ***
## Day:CommunityMultiple               -0.18917    0.08695  -2.176 0.029577 *  
## TreatmentHigh:CommunityMultiple      3.15953    3.23542   0.977 0.328794    
## Day:TreatmentHigh:CommunityMultiple -0.02401    0.12456  -0.193 0.847163    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Day    TrtmnH CmmntM Dy:TrH Dy:CmM TrH:CM
## Day         -0.832                                          
## TreatmntHgh -0.687  0.656                                   
## CmmntyMltpl -0.571  0.525  0.470                            
## Dy:TrtmntHg  0.384 -0.624 -0.838 -0.339                     
## Dy:CmmntyMl  0.494 -0.641 -0.439 -0.918  0.463              
## TrtmntHg:CM  0.390 -0.376 -0.579 -0.742  0.487  0.687       
## Dy:TrtmH:CM -0.278  0.418  0.545  0.624 -0.627 -0.694 -0.928
## convergence code: 0
## Model failed to converge with max|grad| = 0.0634371 (tol = 0.001, component 1)
Anova(survmod3b, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Survivorship
##                           Chisq Df Pr(>Chisq)    
## Day                     60.6800  1  6.715e-15 ***
## Treatment               38.5770  1  5.264e-10 ***
## Community                2.7557  1   0.096909 .  
## Day:Treatment           22.1358  1  2.540e-06 ***
## Day:Community           10.3021  1   0.001329 ** 
## Treatment:Community      4.5938  1   0.032088 *  
## Day:Treatment:Community  0.0371  1   0.847163    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dispersion_glmer(survmod3b) #no overdispersion
## [1] 0.558754
plot(residuals(survmod3b)) + abline(h=0, lty=2) #Dispersed randomly

## integer(0)
meanstresssurv1 <- ddply(final_stress, c("Community", "Treatment"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstresssurv1
##    Community Treatment  N       mean        sd         se max      lower
## 1 Individual   Ambient 21 0.95238095 0.2182179 0.04761905   1  0.7341631
## 2 Individual      High 24 0.00000000 0.0000000 0.00000000   0  0.0000000
## 3   Multiple   Ambient 46 0.76086957 0.4312660 0.06358670   1  0.3296036
## 4   Multiple      High 65 0.04615385 0.2114510 0.02622727   1 -0.1652972
##       upper
## 1 1.1705988
## 2 0.0000000
## 3 1.1921355
## 4 0.2576049

Analyze with full dataset.

#with full time dataset

survmod3<-glmer(Survivorship ~ Day * Treatment * Fusion.Partners * Richness + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family=binomial, subset=c(Community=="Multiple")) 

fixed_form<-glm(Survivorship ~ Day * Treatment * Fusion.Partners, data=stress, family=binomial, subset=c(Community=="Multiple"))

library("brglm2"); glm(fixed_form, data = stress, family = binomial, method="detect_separation")
## Separation: FALSE 
## Existence of maximum likelihood estimates
##                       (Intercept)                               Day 
##                                 0                                 0 
##                     TreatmentHigh                   Fusion.Partners 
##                                 0                                 0 
##                 Day:TreatmentHigh               Day:Fusion.Partners 
##                                 0                                 0 
##     TreatmentHigh:Fusion.Partners Day:TreatmentHigh:Fusion.Partners 
##                                 0                                 0 
## 0: finite value, Inf: infinity, -Inf: -infinity
summary(survmod3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Survivorship ~ Day * Treatment * Fusion.Partners * Richness +  
##     (1 | Parent.Site/Genotype) + (1 | Treatment:Tank) + (1 |      Slide/Sample)
##    Data: stress
##  Subset: c(Community == "Multiple")
## 
##      AIC      BIC   logLik deviance df.resid 
##    266.3    364.0   -112.2    224.3      752 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.4233  0.0000  0.0006  0.0429  1.3089 
## 
## Random effects:
##  Groups               Name        Variance  Std.Dev. 
##  Sample:Slide         (Intercept) 6.397e+00 2.5292012
##  Slide                (Intercept) 1.177e-04 0.0108482
##  Genotype:Parent.Site (Intercept) 8.211e-03 0.0906157
##  Treatment:Tank       (Intercept) 6.207e-04 0.0249143
##  Parent.Site          (Intercept) 1.662e-07 0.0004077
## Number of obs: 773, groups:  
## Sample:Slide, 111; Slide, 45; Genotype:Parent.Site, 20; Treatment:Tank, 8; Parent.Site, 6
## 
## Fixed effects:
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                 9.36537    9.77915   0.958   0.3382
## Day                                        -0.36928    0.33310  -1.109   0.2676
## TreatmentHigh                              -2.26802   10.86374  -0.209   0.8346
## Fusion.Partners                             6.29842    7.40902   0.850   0.3953
## Richness                                    0.91692    7.10758   0.129   0.8974
## Day:TreatmentHigh                           0.02190    0.38514   0.057   0.9546
## Day:Fusion.Partners                        -0.13414    0.24757  -0.542   0.5879
## TreatmentHigh:Fusion.Partners              32.15393   17.20786   1.869   0.0617
## Day:Richness                                0.05145    0.23893   0.215   0.8295
## TreatmentHigh:Richness                      2.42121    7.54725   0.321   0.7484
## Fusion.Partners:Richness                   -1.25492    2.99110  -0.420   0.6748
## Day:TreatmentHigh:Fusion.Partners          -1.42437    0.69415  -2.052   0.0402
## Day:TreatmentHigh:Richness                 -0.16342    0.25928  -0.630   0.5285
## Day:Fusion.Partners:Richness                0.01143    0.10035   0.114   0.9093
## TreatmentHigh:Fusion.Partners:Richness     -5.63393    5.44568  -1.035   0.3009
## Day:TreatmentHigh:Fusion.Partners:Richness  0.26613    0.21337   1.247   0.2123
##                                             
## (Intercept)                                 
## Day                                         
## TreatmentHigh                               
## Fusion.Partners                             
## Richness                                    
## Day:TreatmentHigh                           
## Day:Fusion.Partners                         
## TreatmentHigh:Fusion.Partners              .
## Day:Richness                                
## TreatmentHigh:Richness                      
## Fusion.Partners:Richness                    
## Day:TreatmentHigh:Fusion.Partners          *
## Day:TreatmentHigh:Richness                  
## Day:Fusion.Partners:Richness                
## TreatmentHigh:Fusion.Partners:Richness      
## Day:TreatmentHigh:Fusion.Partners:Richness  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 2 negative eigenvalues
## failure to converge in 10000 evaluations
Anova(survmod3, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Survivorship
##                                          Chisq Df Pr(>Chisq)    
## Day                                    52.1650  1  5.103e-13 ***
## Treatment                              36.7599  1  1.336e-09 ***
## Fusion.Partners                         0.0832  1    0.77299    
## Richness                                1.4063  1    0.23567    
## Day:Treatment                           5.9263  1    0.01492 *  
## Day:Fusion.Partners                     5.5057  1    0.01895 *  
## Treatment:Fusion.Partners               1.0661  1    0.30183    
## Day:Richness                            0.1740  1    0.67658    
## Treatment:Richness                      0.0571  1    0.81121    
## Fusion.Partners:Richness                1.1707  1    0.27927    
## Day:Treatment:Fusion.Partners           5.1440  1    0.02333 *  
## Day:Treatment:Richness                  0.0168  1    0.89683    
## Day:Fusion.Partners:Richness            0.6301  1    0.42733    
## Treatment:Fusion.Partners:Richness      1.3224  1    0.25016    
## Day:Treatment:Fusion.Partners:Richness  1.5557  1    0.21230    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dispersion_glmer(survmod3) #no overdispersion
## [1] 0.4350411
plot(residuals(survmod3)) + abline(h=0, lty=2) #Dispersed randomly

## integer(0)
meanstresssurv2 <- ddply(final_stress, c("Fusion.Partners", "Treatment"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstresssurv2
##   Fusion.Partners Treatment  N       mean        sd         se max      lower
## 1               0   Ambient 30 0.83333333 0.3790490 0.06920457   1  0.4542843
## 2               0      High 50 0.04000000 0.1979487 0.02799417   1 -0.1579487
## 3               1   Ambient 16 0.81250000 0.4031129 0.10077822   1  0.4093871
## 4               1      High 15 0.06666667 0.2581989 0.06666667   1 -0.1915322
## 5               2   Ambient  9 0.77777778 0.4409586 0.14698618   1  0.3368192
## 6               2      High 16 0.00000000 0.0000000 0.00000000   0  0.0000000
## 7               3   Ambient 12 0.83333333 0.3892495 0.11236664   1  0.4440839
## 8               3      High  8 0.00000000 0.0000000 0.00000000   0  0.0000000
##       upper
## 1 1.2123824
## 2 0.2379487
## 3 1.2156129
## 4 0.3248656
## 5 1.2187363
## 6 0.0000000
## 7 1.2225828
## 8 0.0000000

Create plot for effect of number of juveniles (“individual” or “multiple”).

eff.surv3b <- predictorEffect("Community", survmod3b)

surv.effects3b<-plot(eff.surv3b,
                   lines=list(multiline=TRUE, 
                              col=c("blue", "red"), 
                              lty=1), 
                   confint=list(style="bars", width=4, col="black"), 
                   lwd=4,
                   #axes=list(y=list(lim=c(0, 1))),
                   type="rescale", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="none",
                   xlab=expression(bold("Number of Juveniles")), 
                   main="",
                   lattice=list(key.args=list(space="none",
                                              border=FALSE, 
                                              title=expression(bold("Temperature")),
                                              cex=1, 
                                              cex.title=1)));surv.effects3b

Plot effect of number of juveniles over time.

stress$group<-paste(stress$Treatment, "-", stress$Community)

survmod3b2<-glmer(Survivorship ~ group * Day + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide), data=stress, family=binomial, link=logit) 


eff.surv3b2 <- predictorEffect(c("Day"), survmod3b2)
surv.effects3b2<-plot(eff.surv3b2,
                   lines=list(multiline=TRUE, 
                              col=c("blue", "blue", "red", "red"), 
                              lty=c(1, 2, 1, 2)), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1.1))),
                   type="response", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="top",
                   xlab=expression(bold("Day")), 
                   main="",
                   lattice=list(key.args=list(space="top",
                                              border=FALSE, 
                                              title=expression(bold("Temperature - Number of Juveniles")),
                                              cex=1, 
                                              cex.title=1)));surv.effects3b2

Examine curve estimates to identify 50% survivorship.

summary(effect(c("group", "Day"), xlevels=50, survmod3b2))
## 
##  group*Day effect
##                       Day
## group                          0     0.653      1.31      1.96      2.61
##   Ambient - Individual 0.9986417 0.9985680 0.9984898 0.9984082 0.9983223
##   Ambient - Multiple   0.9999603 0.9999528 0.9999437 0.9999331 0.9999204
##   High - Individual    0.9994462 0.9992965 0.9991051 0.9988646 0.9985596
##   High - Multiple      0.9999994 0.9999992 0.9999988 0.9999982 0.9999974
##                       Day
## group                       3.27      3.92      4.57      5.22      5.88
##   Ambient - Individual 0.9982303 0.9981347 0.9980340 0.9979279 0.9978142
##   Ambient - Multiple   0.9999051 0.9998871 0.9998658 0.9998403 0.9998096
##   High - Individual    0.9981660 0.9976737 0.9970497 0.9962588 0.9952396
##   High - Multiple      0.9999961 0.9999943 0.9999917 0.9999879 0.9999823
##                       Day
## group                       6.53      7.18      7.84      8.49      9.14
##   Ambient - Individual 0.9976963 0.9975720 0.9974389 0.9973007 0.9971551
##   Ambient - Multiple   0.9997736 0.9997308 0.9996790 0.9996182 0.9995460
##   High - Individual    0.9939665 0.9923555 0.9902834 0.9877011 0.9844434
##   High - Multiple      0.9999741 0.9999622 0.9999445 0.9999190 0.9998817
##                       Day
## group                        9.8      10.4      11.1      11.8      12.4
##   Ambient - Individual 0.9969992 0.9968502 0.9966669 0.9964730 0.9962979
##   Ambient - Multiple   0.9994586 0.9993647 0.9992345 0.9990775 0.9989175
##   High - Individual    0.9802693 0.9755340 0.9686015 0.9597856 0.9503831
##   High - Multiple      0.9998262 0.9997536 0.9996297 0.9994434 0.9992109
##                       Day
## group                       13.1      13.7      14.4        15      15.7
##   Ambient - Individual 0.9960827 0.9958883 0.9956493 0.9954335 0.9951682
##   Ambient - Multiple   0.9986957 0.9984696 0.9981561 0.9978368 0.9973940
##   High - Individual    0.9367858 0.9224404 0.9019753 0.8807360 0.8510437
##   High - Multiple      0.9988143 0.9983194 0.9974758 0.9964241 0.9946345
##                       Day
## group                       16.3        17      17.6      18.3      18.9
##   Ambient - Individual 0.9949286 0.9946341 0.9943683 0.9940415 0.9937464
##   Ambient - Multiple   0.9969431 0.9963180 0.9956818 0.9948000 0.9939032
##   High - Individual    0.8209593 0.7801006 0.7400643 0.6877657 0.6387043
##   High - Multiple      0.9924079 0.9886311 0.9839535 0.9760725 0.9664033
##                       Day
## group                       19.6      20.2      20.9      21.6      22.2
##   Ambient - Individual 0.9933837 0.9930563 0.9926540 0.9922285 0.9918443
##   Ambient - Multiple   0.9926609 0.9913983 0.9896510 0.9875533 0.9854248
##   High - Individual    0.5776507 0.5232799 0.4592354 0.3965098 0.3452513
##   High - Multiple      0.9503378 0.9310043 0.8997671 0.8565663 0.8080989
##                       Day
## group                       22.9      23.5      24.2      24.8      25.5
##   Ambient - Individual 0.9913724 0.9909463 0.9904229 0.9899504 0.9893700
##   Ambient - Multiple   0.9824856 0.9795085 0.9754059 0.9712606 0.9655648
##   High - Individual    0.2897524 0.2466542 0.2021122 0.1689489 0.1359078
##   High - Multiple      0.7369400 0.6639104 0.5678757 0.4809678 0.3813680
##                       Day
## group                       26.1       26.8       27.4       28.1       28.7
##   Ambient - Individual 0.9888462 0.98820274 0.98762213 0.98690902 0.98626562
##   Ambient - Multiple   0.9598291 0.95197972 0.94411219 0.93340482 0.92274101
##   High - Individual    0.1120816 0.08897129 0.07268148 0.05717206 0.04640781
##   High - Multiple      0.3029899 0.22431753 0.16937894 0.11945350 0.08730690
##                       Day
## group                        29.4         30       30.7       31.3         32
##   Ambient - Individual 0.98547551 0.98476275 0.98388759 0.98309824 0.98212922
##   Ambient - Multiple   0.90833673 0.89411415 0.87509511 0.85652939 0.83202769
##   High - Individual    0.03628552 0.02933141 0.02284452 0.01841715 0.01430846
##   High - Multiple      0.05983017 0.04294646 0.02898720 0.02061633 0.01381049
## 
##  Lower 95 Percent Confidence Limits
##                       Day
## group                          0     0.653      1.31      1.96      2.61
##   Ambient - Individual 0.9646655 0.9649766 0.9652683 0.9655355 0.9657805
##   Ambient - Multiple   0.9979830 0.9977864 0.9975689 0.9973324 0.9970723
##   High - Individual    0.9913895 0.9898742 0.9880798 0.9859913 0.9835369
##   High - Multiple      0.9999769 0.9999692 0.9999588 0.9999452 0.9999270
##                       Day
## group                       3.27      3.92      4.57      5.22      5.88
##   Ambient - Individual 0.9660059 0.9662040 0.9663773 0.9665247 0.9666467
##   Ambient - Multiple   0.9967818 0.9964668 0.9961203 0.9957389 0.9953123
##   High - Individual    0.9806048 0.9772087 0.9732211 0.9685412 0.9629614
##   High - Multiple      0.9999023 0.9998698 0.9998265 0.9997688 0.9996904
##                       Day
## group                       6.53      7.18      7.84      8.49      9.14
##   Ambient - Individual 0.9667381 0.9667996 0.9668295 0.9668250 0.9667847
##   Ambient - Multiple   0.9948491 0.9943388 0.9937671 0.9931458 0.9924601
##   High - Individual    0.9565144 0.9489683 0.9400006 0.9296820 0.9176650
##   High - Multiple      0.9995872 0.9994494 0.9992622 0.9990153 0.9986853
##                       Day
## group                        9.8      10.4      11.1      11.8      12.4
##   Ambient - Individual 0.9667047 0.9665953 0.9664205 0.9661908 0.9659467
##   Ambient - Multiple   0.9916908 0.9909205 0.9899270 0.9888192 0.9877678
##   High - Individual    0.9034731 0.8885987 0.8685655 0.8453005 0.8225148
##   High - Multiple      0.9982364 0.9976957 0.9968508 0.9956941 0.9943680
##                       Day
## group                       13.1      13.7      14.4        15      15.7
##   Ambient - Individual 0.9656022 0.9652515 0.9647720 0.9642954 0.9636556
##   Ambient - Multiple   0.9864083 0.9851157 0.9834407 0.9818446 0.9797712
##   High - Individual    0.7923188 0.7631391 0.7250822 0.6889804 0.6429070
##   High - Multiple      0.9922933 0.9899136 0.9861903 0.9819210 0.9752488
##                       Day
## group                       16.3        17      17.6      18.3      18.9
##   Ambient - Individual 0.9630289 0.9621975 0.9613909 0.9603290 0.9593052
##   Ambient - Multiple   0.9777900 0.9752087 0.9727341 0.9694981 0.9663838
##   High - Individual    0.6002734 0.5474025 0.5000354 0.4433987 0.3946586
##   High - Multiple      0.9676146 0.9557286 0.9422083 0.9213483 0.8979280
##                       Day
## group                       19.6      20.2      20.9      21.6      22.2
##   Ambient - Individual 0.9579646 0.9566780 0.9549996 0.9531103 0.9513053
##   Ambient - Multiple   0.9622932 0.9583378 0.9531159 0.9471326 0.9412965
##   High - Individual    0.3388896 0.2931041 0.2432415 0.1982310 0.1640523
##   High - Multiple      0.8624812 0.8237211 0.7671948 0.6975931 0.6279100
##                       Day
## group                       22.9      23.5       24.2       24.8       25.5
##   Ambient - Individual 0.9489604 0.9467251 0.94382720 0.94107013 0.93750221
##   Ambient - Multiple   0.9335201 0.9258872 0.91565244 0.90554622 0.89192097
##   High - Individual    0.1295765 0.1046011 0.08047226 0.06366912 0.04798639
##   High - Multiple      0.5372850 0.4551979 0.36018275 0.28439130 0.20743679
##                       Day
## group                        26.1       26.8       27.4       28.1       28.7
##   Ambient - Individual 0.93411387 0.92973700 0.92558810 0.92023915 0.91517905
##   Ambient - Multiple   0.87840706 0.86013421 0.84199477 0.81751490 0.79334071
##   High - Individual    0.03739302 0.02775739 0.02139171 0.01570692 0.01200962
##   High - Multiple      0.15343341 0.10461494 0.07369029 0.04796865 0.03274520
##                       Day
## group                         29.4          30        30.7        31.3
##   Ambient - Individual 0.908669500 0.902525769 0.894642221 0.887222149
##   Ambient - Multiple   0.761024773 0.729555825 0.688298487 0.649103361
##   High - Individual    0.008749816 0.006652628 0.004819988 0.003649879
##   High - Multiple      0.020718341 0.013881109 0.008636886 0.005723046
##                       Day
## group                           32
##   Ambient - Individual 0.877729690
##   Ambient - Multiple   0.599268815
##   High - Individual    0.002633742
##   High - Multiple      0.003525362
## 
##  Upper 95 Percent Confidence Limits
##                       Day
## group                          0     0.653      1.31      1.96      2.61
##   Ambient - Individual 0.9999495 0.9999433 0.9999364 0.9999288 0.9999203
##   Ambient - Multiple   0.9999992 0.9999990 0.9999987 0.9999983 0.9999978
##   High - Individual    0.9999647 0.9999516 0.9999335 0.9999091 0.9998757
##   High - Multiple      1.0000000 1.0000000 1.0000000 0.9999999 0.9999999
##                       Day
## group                       3.27      3.92      4.57      5.22      5.88
##   Ambient - Individual 0.9999107 0.9999002 0.9998885 0.9998755 0.9998609
##   Ambient - Multiple   0.9999972 0.9999964 0.9999954 0.9999940 0.9999923
##   High - Individual    0.9998294 0.9997669 0.9996819 0.9995660 0.9994055
##   High - Multiple      0.9999998 0.9999998 0.9999996 0.9999994 0.9999990
##                       Day
## group                       6.53      7.18      7.84      8.49      9.14
##   Ambient - Individual 0.9998451 0.9998275 0.9998079 0.9997866 0.9997631
##   Ambient - Multiple   0.9999901 0.9999873 0.9999836 0.9999789 0.9999728
##   High - Individual    0.9991902 0.9988977 0.9984940 0.9979542 0.9972245
##   High - Multiple      0.9999984 0.9999974 0.9999958 0.9999933 0.9999894
##                       Day
## group                        9.8      10.4      11.1      11.8      12.4
##   Ambient - Individual 0.9997371 0.9997112 0.9996782 0.9996421 0.9996085
##   Ambient - Multiple   0.9999650 0.9999559 0.9999423 0.9999246 0.9999052
##   High - Individual    0.9962224 0.9950079 0.9931037 0.9904985 0.9875264
##   High - Multiple      0.9999829 0.9999737 0.9999566 0.9999283 0.9998899
##                       Day
## group                       13.1      13.7      14.4        15      15.7
##   Ambient - Individual 0.9995660 0.9995267 0.9994773 0.9994319 0.9993753
##   Ambient - Multiple   0.9998762 0.9998445 0.9997974 0.9997459 0.9996695
##   High - Individual    0.9829246 0.9777298 0.9697905 0.9609650 0.9477285
##   High - Multiple      0.9998186 0.9997219 0.9995429 0.9993010 0.9988547
##                       Day
## group                       16.3        17      17.6      18.3      18.9
##   Ambient - Individual 0.9993237 0.9992598 0.9992019 0.9991310 0.9990673
##   Ambient - Multiple   0.9995863 0.9994630 0.9993294 0.9991323 0.9989194
##   High - Individual    0.9333370 0.9123221 0.8901685 0.8589709 0.8273928
##   High - Multiple      0.9982544 0.9971533 0.9956827 0.9930096 0.9894800
##                       Day
## group                       19.6      20.2      20.9      21.6      22.2
##   Ambient - Individual 0.9989901 0.9989215 0.9988391 0.9987546 0.9986808
##   Ambient - Multiple   0.9986070 0.9982714 0.9977819 0.9971622 0.9965044
##   High - Individual    0.7849112 0.7439744 0.6917146 0.6358343 0.5862351
##   High - Multiple      0.9831613 0.9749784 0.9607114 0.9392466 0.9131048
##                       Day
## group                       22.9      23.5      24.2      24.8      25.5
##   Ambient - Individual 0.9985938 0.9985188 0.9984314 0.9983570 0.9982714
##   Ambient - Multiple   0.9955574 0.9945621 0.9931457 0.9916758 0.9896129
##   High - Individual    0.5278548 0.4785237 0.4230312 0.3780295 0.3292136
##   High - Multiple      0.8711125 0.8236434 0.7541648 0.6836186 0.5921726
##                       Day
## group                       26.1      26.8      27.4      28.1      28.7
##   Ambient - Individual 0.9981994 0.9981177 0.9980500 0.9979741 0.9979120
##   Ambient - Multiple   0.9875043 0.9845934 0.9816688 0.9777043 0.9737934
##   High - Individual    0.2908738 0.2504116 0.2193783 0.1872750 0.1630685
##   High - Multiple      0.5104289 0.4171719 0.3432751 0.2675319 0.2127824
##                       Day
## group                       29.4        30       30.7       31.3         32
##   Ambient - Individual 0.9978434 0.9977881 0.99772791 0.99768010 0.99762885
##   Ambient - Multiple   0.9685890 0.9635462 0.95694978 0.95065978 0.94255216
##   High - Individual    0.1383791 0.1199838 0.10140441 0.08767469 0.07389981
##   High - Multiple      0.1606633 0.1251475 0.09279867 0.07148059 0.05252057

Create plot for effect of temperature on survival in multiple juvenile groups across number of fusion partners.

stress$group<-paste(stress$Treatment, "-", stress$Fusion.Partners)

survmod3a<-glmer(Survivorship ~ Day * group + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide), data=stress, family=binomial, link=logit)

summary(effect(c("Day", "group"), xlevels=50, survmod3a))
## 
##  Day*group effect
##        group
## Day     Ambient - 0 Ambient - 1 Ambient - 2 Ambient - 3   High - 0     High - 1
##   0       0.9996572   0.9999989   0.9999970   0.9999981 0.99982061 1.000000e+00
##   0.653   0.9996213   0.9999986   0.9999963   0.9999976 0.99977069 1.000000e+00
##   1.31    0.9995814   0.9999983   0.9999954   0.9999971 0.99970645 1.000000e+00
##   1.96    0.9995378   0.9999978   0.9999942   0.9999964 0.99962520 1.000000e+00
##   2.61    0.9994896   0.9999972   0.9999928   0.9999955 0.99952147 1.000000e+00
##   3.27    0.9994355   0.9999964   0.9999911   0.9999944 0.99938676 1.000000e+00
##   3.92    0.9993767   0.9999955   0.9999889   0.9999931 0.99921710 1.000000e+00
##   4.57    0.9993118   0.9999943   0.9999862   0.9999914 0.99900055 1.000000e+00
##   5.22    0.9992400   0.9999927   0.9999828   0.9999894 0.99872418 1.000000e+00
##   5.88    0.9991596   0.9999907   0.9999785   0.9999868 0.99836538 1.000000e+00
##   6.53    0.9990720   0.9999883   0.9999733   0.9999837 0.99791373 1.000000e+00
##   7.18    0.9989753   0.9999851   0.9999668   0.9999798 0.99733763 1.000000e+00
##   7.84    0.9988669   0.9999811   0.9999585   0.9999749 0.99659022 1.000000e+00
##   8.49    0.9987489   0.9999760   0.9999484   0.9999690 0.99565023 1.000000e+00
##   9.14    0.9986186   0.9999696   0.9999359   0.9999616 0.99445256 1.000000e+00
##   9.8     0.9984724   0.9999614   0.9999200   0.9999522 0.99290099 9.999999e-01
##   10.4    0.9983262   0.9999519   0.9999021   0.9999418 0.99112022 9.999999e-01
##   11.1    0.9981378   0.9999379   0.9998762   0.9999268 0.98847721 9.999997e-01
##   11.8    0.9979283   0.9999199   0.9998435   0.9999078 0.98505939 9.999993e-01
##   12.4    0.9977301   0.9999003   0.9998086   0.9998877 0.98134852 9.999985e-01
##   13.1    0.9974748   0.9998713   0.9997579   0.9998587 0.97586784 9.999964e-01
##   13.7    0.9972334   0.9998399   0.9997040   0.9998278 0.96994367 9.999923e-01
##   14.4    0.9969224   0.9997934   0.9996256   0.9997833 0.96124408 9.999818e-01
##   15      0.9966283   0.9997428   0.9995422   0.9997360 0.95190686 9.999616e-01
##   15.7    0.9962496   0.9996681   0.9994210   0.9996677 0.93831863 9.999087e-01
##   16.3    0.9958914   0.9995870   0.9992920   0.9995953 0.92389542 9.998079e-01
##   17      0.9954303   0.9994670   0.9991047   0.9994906 0.90319768 9.995427e-01
##   17.6    0.9949943   0.9993367   0.9989053   0.9993796 0.88159829 9.990384e-01
##   18.3    0.9944331   0.9991441   0.9986159   0.9992192 0.85124931 9.977130e-01
##   18.9    0.9939025   0.9989350   0.9983077   0.9990490 0.82036409 9.952009e-01
##   19.6    0.9932198   0.9986258   0.9978607   0.9988032 0.77826709 9.886469e-01
##   20.2    0.9925744   0.9982904   0.9973848   0.9985426 0.73691184 9.764115e-01
##   20.9    0.9917442   0.9977945   0.9966948   0.9981661 0.68281940 9.456000e-01
##   21.6    0.9908219   0.9971551   0.9958234   0.9976925 0.62329066 8.795086e-01
##   22.2    0.9899507   0.9964619   0.9948967   0.9971907 0.56903668 7.762713e-01
##   22.9    0.9888303   0.9954380   0.9935543   0.9964661 0.50367541 5.930040e-01
##   23.5    0.9877724   0.9943289   0.9921282   0.9956987 0.44746580 4.091888e-01
##   24.2    0.9864125   0.9926922   0.9900649   0.9945915 0.38363796 2.253093e-01
##   24.8    0.9851290   0.9909216   0.9878761   0.9934198 0.33186713 1.214567e-01
##   25.5    0.9834801   0.9883132   0.9847154   0.9917308 0.27628347 5.486888e-02
##   26.1    0.9819246   0.9854971   0.9813702   0.9899455 0.23351129 2.685464e-02
##   26.8    0.9799276   0.9813593   0.9765538   0.9873761 0.18972323 1.145549e-02
##   27.4    0.9780451   0.9769067   0.9714741   0.9846651 0.15743652 5.478218e-03
##   28.1    0.9756302   0.9703918   0.9641925   0.9807725 0.12557683 2.307795e-03
##   28.7    0.9733556   0.9634171   0.9565539   0.9766769 0.10282118 1.098325e-03
##   29.4    0.9704405   0.9532786   0.9456767   0.9708167 0.08095186 4.615125e-04
##   30      0.9676976   0.9425114   0.9343572   0.9646768 0.06567531 2.194302e-04
##   30.7    0.9641861   0.9270180   0.9183968   0.9559376 0.05125529 9.215677e-05
##   31.3    0.9608862   0.9107653   0.9019821   0.9468387 0.04133080 4.380825e-05
##   32      0.9566675   0.8877354   0.8791680   0.9339881 0.03207250 1.839681e-05
##        group
## Day         High - 2     High - 3
##   0     1.000000e+00 1.000000e+00
##   0.653 1.000000e+00 1.000000e+00
##   1.31  1.000000e+00 1.000000e+00
##   1.96  1.000000e+00 1.000000e+00
##   2.61  1.000000e+00 1.000000e+00
##   3.27  1.000000e+00 1.000000e+00
##   3.92  1.000000e+00 1.000000e+00
##   4.57  1.000000e+00 1.000000e+00
##   5.22  1.000000e+00 1.000000e+00
##   5.88  1.000000e+00 1.000000e+00
##   6.53  1.000000e+00 1.000000e+00
##   7.18  1.000000e+00 1.000000e+00
##   7.84  1.000000e+00 1.000000e+00
##   8.49  1.000000e+00 1.000000e+00
##   9.14  1.000000e+00 1.000000e+00
##   9.8   1.000000e+00 1.000000e+00
##   10.4  1.000000e+00 1.000000e+00
##   11.1  1.000000e+00 1.000000e+00
##   11.8  1.000000e+00 1.000000e+00
##   12.4  1.000000e+00 1.000000e+00
##   13.1  1.000000e+00 1.000000e+00
##   13.7  1.000000e+00 1.000000e+00
##   14.4  1.000000e+00 9.999999e-01
##   15    1.000000e+00 9.999999e-01
##   15.7  1.000000e+00 9.999996e-01
##   16.3  1.000000e+00 9.999989e-01
##   17    9.999999e-01 9.999967e-01
##   17.6  9.999996e-01 9.999915e-01
##   18.3  9.999982e-01 9.999744e-01
##   18.9  9.999937e-01 9.999345e-01
##   19.6  9.999719e-01 9.998034e-01
##   20.2  9.998995e-01 9.994961e-01
##   20.9  9.995553e-01 9.984897e-01
##   21.6  9.980336e-01 9.954827e-01
##   22.2  9.929953e-01 9.884978e-01
##   22.9  9.697086e-01 9.662694e-01
##   23.5  8.994109e-01 9.178410e-01
##   24.2  6.687835e-01 7.883074e-01
##   24.8  3.606032e-01 5.922039e-01
##   25.5  1.129702e-01 3.261762e-01
##   26.1  3.435032e-02 1.587978e-01
##   26.8  7.968983e-03 5.919969e-02
##   27.4  2.238669e-03 2.395144e-02
##   28.1  5.064186e-04 8.113347e-03
##   28.7  1.414991e-04 3.179750e-03
##   29.4  3.195713e-05 1.062166e-03
##   30    8.926134e-06 4.144881e-04
##   30.7  2.015732e-06 1.382006e-04
##   31.3  5.630138e-07 5.389950e-05
##   32    1.271410e-07 1.796711e-05
## 
##  Lower 95 Percent Confidence Limits
##        group
## Day     Ambient - 0 Ambient - 1 Ambient - 2 Ambient - 3    High - 0
##   0       0.9912518   0.9545168   0.7325815   0.8462365 0.997701705
##   0.653   0.9908699   0.9539437   0.7364089   0.8484868 0.997242915
##   1.31    0.9904657   0.9533518   0.7401424   0.8506719 0.996688404
##   1.96    0.9900450   0.9527502   0.7437172   0.8527544 0.996029605
##   2.61    0.9896022   0.9521316   0.7471703   0.8547566 0.995239134
##   3.27    0.9891284   0.9514850   0.7505482   0.8567057 0.994274598
##   3.92    0.9886362   0.9508288   0.7537444   0.8585408 0.993133072
##   4.57    0.9881170   0.9501518   0.7568066   0.8602898 0.991763107
##   5.22    0.9875687   0.9494527   0.7597298   0.8619502 0.990118988
##   5.88    0.9869801   0.9487184   0.7625502   0.8635426 0.988112714
##   6.53    0.9863666   0.9479692   0.7651760   0.8650155 0.985738645
##   7.18    0.9857169   0.9471922   0.7676444   0.8663899 0.982890630
##   7.84    0.9850174   0.9463723   0.7699816   0.8676802 0.979417641
##   8.49    0.9842859   0.9455317   0.7721083   0.8688425 0.975312202
##   9.14    0.9835088   0.9446553   0.7740511   0.8698915 0.970394334
##   9.8     0.9826691   0.9437253   0.7758245   0.8708339 0.964409325
##   10.4    0.9818578   0.9428414   0.7772512   0.8715769 0.957941763
##   11.1    0.9808483   0.9417593   0.7786770   0.8722975 0.948930606
##   11.8    0.9797648   0.9406163   0.7798260   0.8728488 0.938048523
##   12.4    0.9787716   0.9395827   0.7805725   0.8731755 0.926966414
##   13.1    0.9775303   0.9383067   0.7811410   0.8733718 0.911640875
##   13.7    0.9763890   0.9371458   0.7813462   0.8733674 0.896137579
##   14.4    0.9749581   0.9357031   0.7812235   0.8731403 0.874877034
##   15      0.9736381   0.9343812   0.7807770   0.8727359 0.853587459
##   15.7    0.9719774   0.9327256   0.7798132   0.8719913 0.824756426
##   16.3    0.9704400   0.9311960   0.7785645   0.8710918 0.796316308
##   17      0.9684988   0.9292627   0.7765530   0.8696980 0.758497888
##   17.6    0.9666950   0.9274589   0.7742936   0.8681694 0.721983661
##   18.3    0.9644087   0.9251547   0.7709464   0.8659395 0.674651080
##   18.9    0.9622761   0.9229797   0.7673835   0.8635894 0.630273632
##   19.6    0.9595626   0.9201658   0.7622936   0.8602527 0.574670885
##   20.2    0.9570216   0.9174732   0.7570106   0.8568006 0.524491444
##   20.9    0.9537761   0.9139373   0.7495969   0.8519589 0.464245411
##   21.6    0.9501908   0.9098839   0.7406070   0.8460764 0.403748426
##   22.2    0.9468142   0.9059037   0.7314288   0.8400456 0.353041394
##   22.9    0.9424772   0.9005273   0.7186987   0.8316244 0.296738636
##   23.5    0.9383789   0.8951399   0.7057623   0.8229870 0.251991951
##   24.2    0.9330984   0.8877021   0.6878924   0.8109048 0.204863053
##   24.8    0.9280945   0.8800777   0.6698116   0.7984833 0.169291314
##   25.5    0.9216313   0.8692986   0.6449788   0.7810695 0.133604370
##   26.1    0.9154940   0.8579844   0.6200568   0.7631454 0.107867452
##   26.8    0.9075543   0.8416163   0.5862431   0.7380417 0.083081857
##   27.4    0.9000072   0.8240793   0.5528906   0.7123116 0.065853898
##   28.1    0.8902394   0.7982953   0.5087392   0.6766100 0.049782817
##   28.7    0.8809566   0.7704196   0.4665948   0.6405886 0.038919461
##   29.4    0.8689542   0.7295108   0.4131631   0.5918255 0.029020218
##   30      0.8575680   0.6859797   0.3648143   0.5442980 0.022462141
##   30.7    0.8428857   0.6242950   0.3073540   0.4829192 0.016584592
##   31.3    0.8290071   0.5622163   0.2590662   0.4265387 0.012745536
##   32      0.8111942   0.4810510   0.2061730   0.3588126 0.009344534
##        group
## Day         High - 1     High - 2     High - 3
##   0     9.999827e-01 1.000000e+00 9.998425e-01
##   0.653 9.999758e-01 1.000000e+00 9.998034e-01
##   1.31  9.999662e-01 1.000000e+00 9.997543e-01
##   1.96  9.999530e-01 1.000000e+00 9.996936e-01
##   2.61  9.999346e-01 1.000000e+00 9.996177e-01
##   3.27  9.999084e-01 1.000000e+00 9.995213e-01
##   3.92  9.998725e-01 1.000000e+00 9.994025e-01
##   4.57  9.998223e-01 1.000000e+00 9.992540e-01
##   5.22  9.997525e-01 1.000000e+00 9.990682e-01
##   5.88  9.996533e-01 9.999999e-01 9.988318e-01
##   6.53  9.995167e-01 9.999999e-01 9.985398e-01
##   7.18  9.993260e-01 9.999998e-01 9.981740e-01
##   7.84  9.990552e-01 9.999996e-01 9.977077e-01
##   8.49  9.986818e-01 9.999993e-01 9.971307e-01
##   9.14  9.981604e-01 9.999987e-01 9.964066e-01
##   9.8   9.974187e-01 9.999977e-01 9.954813e-01
##   10.4  9.964869e-01 9.999960e-01 9.944316e-01
##   11.1  9.949653e-01 9.999924e-01 9.928895e-01
##   11.8  9.927832e-01 9.999856e-01 9.909120e-01
##   12.4  9.901732e-01 9.999750e-01 9.887758e-01
##   13.1  9.859144e-01 9.999524e-01 9.856261e-01
##   13.7  9.808288e-01 9.999172e-01 9.822144e-01
##   14.4  9.725550e-01 9.998419e-01 9.771691e-01
##   15    9.627213e-01 9.997243e-01 9.716861e-01
##   15.7  9.468436e-01 9.994717e-01 9.635476e-01
##   16.3  9.281818e-01 9.990757e-01 9.546667e-01
##   17    8.985543e-01 9.982206e-01 9.414220e-01
##   17.6  8.645462e-01 9.968729e-01 9.268921e-01
##   18.3  8.123640e-01 9.939432e-01 9.050880e-01
##   18.9  7.551412e-01 9.892942e-01 8.809994e-01
##   19.6  6.726272e-01 9.791248e-01 8.445556e-01
##   20.2  5.889855e-01 9.629366e-01 8.039375e-01
##   20.9  4.797506e-01 9.277187e-01 7.419533e-01
##   21.6  3.648355e-01 8.606862e-01 6.591930e-01
##   22.2  2.690919e-01 7.616605e-01 5.673259e-01
##   22.9  1.698139e-01 5.817398e-01 4.331900e-01
##   23.5  1.020533e-01 3.838546e-01 3.015108e-01
##   24.2  4.799674e-02 1.695646e-01 1.563212e-01
##   24.8  2.183644e-02 5.899768e-02 6.913025e-02
##   25.5  7.583508e-03 1.161293e-02 1.995638e-02
##   26.1  2.808945e-03 2.266247e-03 5.660537e-03
##   26.8  8.250626e-04 2.859700e-04 1.117317e-03
##   27.4  2.780836e-04 4.474087e-05 2.552580e-04
##   28.1  7.598612e-05 4.866239e-06 4.274697e-05
##   28.7  2.457363e-05 7.058022e-07 8.901133e-06
##   29.4  6.494592e-06 7.260488e-08 1.384993e-06
##   30    2.058156e-06 1.020469e-08 2.760031e-07
##   30.7  5.346023e-07 1.023491e-09 4.139013e-08
##   31.3  1.675413e-07 1.416432e-10 8.058400e-09
##   32    4.308669e-08 1.401818e-11 1.183879e-09
## 
##  Upper 95 Percent Confidence Limits
##        group
## Day     Ambient - 0 Ambient - 1 Ambient - 2 Ambient - 3  High - 0    High - 1
##   0       0.9999867   1.0000000   1.0000000   1.0000000 0.9999860 1.000000000
##   0.653   0.9999844   1.0000000   1.0000000   1.0000000 0.9999810 1.000000000
##   1.31    0.9999818   1.0000000   1.0000000   1.0000000 0.9999740 1.000000000
##   1.96    0.9999787   1.0000000   1.0000000   1.0000000 0.9999647 1.000000000
##   2.61    0.9999752   1.0000000   1.0000000   1.0000000 0.9999521 1.000000000
##   3.27    0.9999710   1.0000000   1.0000000   1.0000000 0.9999346 1.000000000
##   3.92    0.9999662   1.0000000   1.0000000   1.0000000 0.9999112 1.000000000
##   4.57    0.9999606   1.0000000   1.0000000   1.0000000 0.9998795 1.000000000
##   5.22    0.9999541   1.0000000   1.0000000   1.0000000 0.9998365 1.000000000
##   5.88    0.9999464   1.0000000   1.0000000   1.0000000 0.9997772 1.000000000
##   6.53    0.9999376   1.0000000   1.0000000   1.0000000 0.9996980 1.000000000
##   7.18    0.9999274   1.0000000   1.0000000   1.0000000 0.9995908 1.000000000
##   7.84    0.9999154   1.0000000   1.0000000   1.0000000 0.9994433 1.000000000
##   8.49    0.9999017   1.0000000   1.0000000   1.0000000 0.9992466 1.000000000
##   9.14    0.9998859   1.0000000   1.0000000   1.0000000 0.9989811 1.000000000
##   9.8     0.9998673   1.0000000   1.0000000   1.0000000 0.9986167 1.000000000
##   10.4    0.9998479   1.0000000   1.0000000   1.0000000 0.9981751 1.000000000
##   11.1    0.9998218   0.9999999   0.9999999   1.0000000 0.9974814 1.000000000
##   11.8    0.9997914   0.9999999   0.9999999   0.9999999 0.9965288 1.000000000
##   12.4    0.9997614   0.9999998   0.9999999   0.9999999 0.9954361 1.000000000
##   13.1    0.9997213   0.9999997   0.9999998   0.9999999 0.9937302 0.999999999
##   13.7    0.9996818   0.9999996   0.9999997   0.9999998 0.9917830 0.999999997
##   14.4    0.9996291   0.9999994   0.9999995   0.9999997 0.9887614 0.999999988
##   15      0.9995775   0.9999991   0.9999993   0.9999995 0.9853367 0.999999962
##   15.7    0.9995087   0.9999985   0.9999988   0.9999992 0.9800682 0.999999851
##   16.3    0.9994416   0.9999977   0.9999982   0.9999989 0.9741575 0.999999523
##   17      0.9993525   0.9999963   0.9999972   0.9999983 0.9651785 0.999998146
##   17.6    0.9992659   0.9999944   0.9999959   0.9999975 0.9552544 0.999994086
##   18.3    0.9991516   0.9999909   0.9999935   0.9999961 0.9404515 0.999977252
##   18.9    0.9990409   0.9999864   0.9999905   0.9999943 0.9244386 0.999928290
##   19.6    0.9988954   0.9999782   0.9999853   0.9999912 0.9011669 0.999729128
##   20.2    0.9987553   0.9999674   0.9999786   0.9999873 0.8767399 0.999164361
##   20.9    0.9985722   0.9999481   0.9999671   0.9999806 0.8424778 0.996957269
##   21.6    0.9983658   0.9999178   0.9999498   0.9999706 0.8016993 0.989334387
##   22.2    0.9981689   0.9998786   0.9999283   0.9999583 0.7616126 0.970326274
##   22.9    0.9979138   0.9998099   0.9998925   0.9999379 0.7093606 0.912115385
##   23.5    0.9976719   0.9997224   0.9998490   0.9999132 0.6606486 0.808450940
##   24.2    0.9973606   0.9995718   0.9997781   0.9998732 0.6005848 0.626552626
##   24.8    0.9970674   0.9993844   0.9996945   0.9998262 0.5476434 0.461248455
##   25.5    0.9966928   0.9990708   0.9995625   0.9997520 0.4858822 0.306063530
##   26.1    0.9963424   0.9986933   0.9994122   0.9996677 0.4342658 0.212812347
##   26.8    0.9958979   0.9980865   0.9991839   0.9995397 0.3769704 0.139877979
##   27.4    0.9954850   0.9973892   0.9989349   0.9993998 0.3312229 0.098353582
##   28.1    0.9949650   0.9963290   0.9985738   0.9991965 0.2824639 0.065778486
##   28.7    0.9944854   0.9951847   0.9981987   0.9989846 0.2449067 0.046889950
##   29.4    0.9938854   0.9935632   0.9976822   0.9986915 0.2060904 0.031782362
##   30      0.9933358   0.9919383   0.9971732   0.9984011 0.1769724 0.022869494
##   30.7    0.9926526   0.9898059   0.9965089   0.9980197 0.1475328 0.015640733
##   31.3    0.9920307   0.9878220   0.9958880   0.9976607 0.1258535 0.011326114
##   32      0.9912622   0.9853919   0.9951180   0.9972124 0.1042620 0.007793989
##        group
## Day        High - 2  High - 3
##   0     1.000000000 1.0000000
##   0.653 1.000000000 1.0000000
##   1.31  1.000000000 1.0000000
##   1.96  1.000000000 1.0000000
##   2.61  1.000000000 1.0000000
##   3.27  1.000000000 1.0000000
##   3.92  1.000000000 1.0000000
##   4.57  1.000000000 1.0000000
##   5.22  1.000000000 1.0000000
##   5.88  1.000000000 1.0000000
##   6.53  1.000000000 1.0000000
##   7.18  1.000000000 1.0000000
##   7.84  1.000000000 1.0000000
##   8.49  1.000000000 1.0000000
##   9.14  1.000000000 1.0000000
##   9.8   1.000000000 1.0000000
##   10.4  1.000000000 1.0000000
##   11.1  1.000000000 1.0000000
##   11.8  1.000000000 1.0000000
##   12.4  1.000000000 1.0000000
##   13.1  1.000000000 1.0000000
##   13.7  1.000000000 1.0000000
##   14.4  1.000000000 1.0000000
##   15    1.000000000 1.0000000
##   15.7  1.000000000 1.0000000
##   16.3  1.000000000 1.0000000
##   17    1.000000000 1.0000000
##   17.6  1.000000000 1.0000000
##   18.3  0.999999999 1.0000000
##   18.9  0.999999996 1.0000000
##   19.6  0.999999963 0.9999998
##   20.2  0.999999738 0.9999990
##   20.9  0.999997459 0.9999934
##   21.6  0.999976017 0.9999602
##   22.2  0.999841004 0.9998225
##   22.9  0.998644654 0.9990696
##   23.5  0.992267885 0.9965532
##   24.2  0.952306511 0.9868145
##   24.8  0.835339624 0.9659833
##   25.5  0.579920206 0.9200471
##   26.1  0.357778717 0.8622565
##   26.8  0.184063272 0.7797273
##   27.4  0.101134071 0.7022452
##   28.1  0.050111396 0.6101571
##   28.7  0.027592771 0.5333971
##   29.4  0.013871742 0.4494359
##   30    0.007747416 0.3838474
##   30.7  0.003954239 0.3158071
##   31.3  0.002232915 0.2650040
##   32    0.001151806 0.2142610
eff.surv3 <- predictorEffect(c("Day"), survmod3a)
surv.effects3<-plot(eff.surv3,
                   lines=list(multiline=TRUE, 
                              col=c("blue", "blue", "blue", "blue", "red", "red", "red", "red"), 
                              lty=c(1, 2, 4, 3, 1, 2, 4, 3)), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1.1))),
                   type="response", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="top",
                   xlab=expression(bold("Day")), 
                   main="",
                   lattice=list(key.args=list(space="top",
                                              border=FALSE, 
                                              title=expression(bold("Temperature - Fusion Partners")),
                                              cex=1, 
                                              cex.title=1)));surv.effects3

Display mean survivorship within multiple juvenile groups at the end of the experiment.

meanstresssurv2 <- ddply(multiple_final_stress, c("Fusion.Partners", "Treatment"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstresssurv2
##   Fusion.Partners Treatment  N       mean        sd         se max       lower
## 1               0   Ambient  9 0.55555556 0.5270463 0.17568209   1  0.02850928
## 2               0      High 26 0.07692308 0.2717465 0.05329387   1 -0.19482341
## 3               1   Ambient 16 0.81250000 0.4031129 0.10077822   1  0.40938711
## 4               1      High 15 0.06666667 0.2581989 0.06666667   1 -0.19153222
## 5               2   Ambient  9 0.77777778 0.4409586 0.14698618   1  0.33681923
## 6               2      High 16 0.00000000 0.0000000 0.00000000   0  0.00000000
## 7               3   Ambient 12 0.83333333 0.3892495 0.11236664   1  0.44408386
## 8               3      High  8 0.00000000 0.0000000 0.00000000   0  0.00000000
##       upper
## 1 1.0826018
## 2 0.3486696
## 3 1.2156129
## 4 0.3248656
## 5 1.2187363
## 6 0.0000000
## 7 1.2225828
## 8 0.0000000

5. STRESS FUSION

Analysis

Analyze with a logistic binomial regression model.

Analyze fusion success over the course of the stress test.

#with full
fusemod2<-glmer(Fusion ~  Day * Treatment * Richness + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family=binomial, link=logit, subset=c(Community=="Multiple"))

summary(fusemod2)
Anova(fusemod2, type="II") 

plot(allEffects(fusemod2))

dispersion_glmer(fusemod2) #no overdispersion
plot(residuals(fusemod2)) + abline(h=0, lty=2)

There is an effect of treatment and day and richness. Plot this relationship.

stress$group<-paste(stress$Treatment, "-", stress$Richness)

fusemod2a<-glmer(Fusion ~  Day * group + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family=binomial, link=logit, subset=c(Community=="Multiple"))

eff.fuse2 <- predictorEffect("Day", xlevels=4, fusemod2a, rug=TRUE)

fuse.effects2<-plot(eff.fuse2,
                   lines=list(multiline=TRUE, 
                              col=c("blue", "blue", "blue", "blue", "red", "red", "red", "red"), 
                              lty=c(1, 2, 4, 3, 1, 2, 4, 3)),
                   confint=list(style="none"), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Successful Fusion)")), 
                   legend.position="right",
                   xlab=expression(bold("Time (Days)")), 
                   main="",
                   lattice=list(key.args=list(space="none",
                                              border=FALSE, 
                                              title=expression(bold("Temperature - \nGenotypic Richness")),
                                              cex=1, 
                                              cex.title=1)));fuse.effects2

fusemod2b<-glmer(Fusion ~  Day * Treatment + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family=binomial, link=logit, subset=c(Community=="Multiple"))
Anova(fusemod2b, type="II") 
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Fusion
##                 Chisq Df Pr(>Chisq)    
## Day           32.6383  1   1.11e-08 ***
## Treatment      0.2578  1     0.6117    
## Day:Treatment  0.0725  1     0.7877    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
eff.fuse2b <- predictorEffect("Day", xlevels=4, fusemod2b, rug=TRUE)

fuse.effects2b<-plot(eff.fuse2b,
                   lines=list(multiline=TRUE, 
                              col=c("blue", "red"), 
                              lty=c(1)),
                   confint=list(style="bands"), 
                   lwd=4,
                   axes=list(y=list(lim=c(0.8, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Successful Fusion)")), 
                   legend.position="right",
                   xlab=expression(bold("Time (Days)")), 
                   main="",
                   lattice=list(key.args=list(space="none",
                                              border=FALSE, 
                                              title=expression(bold("Temperature")),
                                              cex=1, 
                                              cex.title=1)));fuse.effects2b

Analyze fusion as a product of treatment and genotypic richness at the end of the experiment.

#with final dataset
fusemod3<-glmer(Fusion ~ Treatment   + (1|Treatment:Tank) + (1|Parent.Site/Genotype) + (1|Slide), data=stress, link=logit, family=binomial, subset=c(Community=="Multiple" & Timepoint == "End"))
summary(fusemod3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Fusion ~ Treatment + (1 | Treatment:Tank) + (1 | Parent.Site/Genotype) +  
##     (1 | Slide)
##    Data: stress
##  Subset: c(Community == "Multiple" & Timepoint == "End")
## 
##      AIC      BIC   logLik deviance df.resid 
##    117.0    133.2    -52.5    105.0      105 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.43276 -0.20590  0.04881  0.15110  1.32248 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  Slide                (Intercept) 40.4085  6.3568  
##  Genotype:Parent.Site (Intercept)  0.8975  0.9474  
##  Treatment:Tank       (Intercept)  0.0000  0.0000  
##  Parent.Site          (Intercept)  0.0000  0.0000  
## Number of obs: 111, groups:  
## Slide, 45; Genotype:Parent.Site, 20; Treatment:Tank, 8; Parent.Site, 6
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept)      5.532      2.268   2.439   0.0147 *
## TreatmentHigh   -4.525      2.568  -1.762   0.0781 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TreatmntHgh -0.658
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(fusemod3, type="II") # no effects on rate of fusion by the end point
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Fusion
##            Chisq Df Pr(>Chisq)  
## Treatment 3.1045  1    0.07808 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(fusemod3))

dispersion_glmer(fusemod3) #no overdispersion
## [1] 0.6023656
plot(residuals(fusemod3)) + abline(h=0, lty=2)

## integer(0)
meanstressfuse <- ddply(multiple_final_stress, c("Treatment", "Richness", "Timepoint"), summarise,
                   N    = length(Fusion[!is.na(Fusion)]),
                   mean = mean(Fusion, na.rm=TRUE),
                   sd   = sd(Fusion, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Fusion, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstressfuse
##   Treatment Richness Timepoint  N      mean        sd        se max       lower
## 1   Ambient        1       End 17 0.7058824 0.4696682 0.1139113   1  0.23621413
## 2   Ambient        2       End 14 0.7142857 0.4688072 0.1252940   1  0.24547848
## 3   Ambient        3       End  3 1.0000000 0.0000000 0.0000000   1  1.00000000
## 4   Ambient        4       End 12 1.0000000 0.0000000 0.0000000   1  1.00000000
## 5      High        1       End 25 0.6000000 0.5000000 0.1000000   1  0.10000000
## 6      High        2       End  6 0.6666667 0.5163978 0.2108185   1  0.15026889
## 7      High        3       End 18 0.5000000 0.5144958 0.1212678   1 -0.01449576
## 8      High        4       End 16 0.6875000 0.4787136 0.1196784   1  0.20878645
##      upper
## 1 1.175551
## 2 1.183093
## 3 1.000000
## 4 1.000000
## 5 1.100000
## 6 1.183064
## 7 1.014496
## 8 1.166214
eff.fuse3 <- predictorEffect("Treatment", fusemod3, xlevels=4, rug=TRUE)

fuse.effects3<-plot(eff.fuse3,
                   lines=list(multiline=TRUE, 
                              col=c("black"), 
                              lty=c(1,2,4,3)),
                   confint=list(style="bars"), 
                   lwd=4,
                   #axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Successful Fusion)")), 
                   legend.position="right",
                   xlab=expression(bold("Temperature Treatment")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Genotypic \nRichness")),
                                              cex=1, 
                                              cex.title=1)));fuse.effects3

No effects at the end of the experiment, we will report full analysis over time during the stress test.

6. STRESS GROWTH

Analyze with a poisson regression model.

Analyze growth during the course of the stress test.

Test for effect of temperature:

#with full dataset 
#polypmod2<-glmer(Polyps ~ Treatment + Community + Community:Treatment + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide), data=stress, family = poisson, subset=c(Timepoint=="S2D11"))

#summary(polypmod2)
#Anova(polypmod2, type="II")
#qqPlot(residuals(polypmod2))

polypmod2a<-glmer(Polyps^2 ~ Day * Treatment * Community + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family = poisson)

summary(polypmod2a)
Anova(polypmod2a, type="II")
qqPlot(residuals(polypmod2a))

Plot difference in temperature:

stress$group<-paste(stress$Treatment, "-", stress$Community)

polypmod2b<-glmer(Polyps^2 ~ Day * group + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family = poisson)

eff.surv5 <- predictorEffect("Day", xlevels=4, polypmod2b)

polyp.effects2<-plot(eff.surv5,
                   lines=list(multiline=TRUE, 
                              col=c("blue", "blue", "red", "red"), 
                              lty=c(1,2, 1, 2)), 
                   confint=list(style="bands", width=4), 
                   lwd=4,
                   #axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("(Total Polyp Expansion)"^2)), 
                   legend.position="right",
                   xlab=expression(bold("Time (Days)")), 
                   main="", 
                   lattice=list(key.args=list(space="top",
                                              border=FALSE, 
                                              title=expression(bold("Treatment - Number of Juveniles")),
                                              cex=1, 
                                              cex.title=1)));polyp.effects2

Present means of growth by number of juveniles. Sample size was too low due to mortality at the end of the experiment to calculate growth within multiple juvenile groups, so we pool these together for analysis here.

meanstressgrowth <- ddply(stress, c("Community", "Treatment", "Timepoint"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Polyps, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstressgrowth

Combine figures for Stress Period

stress_figs<-plot_grid(fuse.effects2, surv.effects3b2, surv.effects3, polyp.effects2, labels = c("A", "B", "C", "D"), ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(0.6,1,1,1), label_size = 20, label_y=1, align="h");stress_figs

ggsave(filename="Figures/stress_figs.png", plot=stress_figs, dpi=500, width=20, height=6, units="in")

7. TEMPERATURE DATA

Load in temperature data from Apex files. Data were recorded every 15 minutes and probes were calibrated to a digital thermometer weekly.

juvtemps<-read.csv("Data/ApexTemps.csv", header=T, na.strings="NA")
#convert Date and Time to date/time format
juvtemps$Date <- as.POSIXct(juvtemps$Date, format="%m/%d/%y %H:%M")

#convert to long format
juvtemps<-juvtemps %>% gather(Tank, Temperature, Tank17:Tank24)
juvtemps$Tank<-as.factor(juvtemps$Tank)

dailymax<-juvtemps[which(juvtemps$TimeDay == "930"),]
dailymax<-juvtemps[which(juvtemps$Phase == "Stress"),]

dailymax$Treatment<-ifelse(dailymax$Tank %in% c("Tank17"), "Ambient",
                            ifelse(dailymax$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(dailymax$Tank %in% c("Tank19"), "High",
                            ifelse(dailymax$Tank %in% c("Tank20"), "High",
                            ifelse(dailymax$Tank %in% c("Tank21"), "Ambient",
                            ifelse(dailymax$Tank %in% c("Tank22"), "High",
                            ifelse(dailymax$Tank %in% c("Tank23"), "Ambient",
                            ifelse(dailymax$Tank %in% c("Tank24"), "High",NA))))))))

maxtemps <- ddply(dailymax, c("Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Temperature, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )

juvtemps$Period<-ifelse(juvtemps$Phase %in% c("Stress"), "Stress",
                            ifelse(juvtemps$Phase %in% c("Stress2"), "Stress", 
                            ifelse(juvtemps$Phase %in% c("Acclimation"), "Stress",
                            ifelse(juvtemps$Phase %in% c("Hurricane"), "Hurricane",
                            ifelse(juvtemps$Phase %in% c("Ambient"), "Grow-Out",NA)))))

Plot timeseries of temperature data with tank temperature in colors.

#subset temperature measurements

juvtempsHIGH<-subset(juvtemps, Period==c("Stress"), c(Date, TimeDay, Phase, Tank, Temperature, Period))
juvtempsHIGH$Treatment<-ifelse(juvtempsHIGH$Tank %in% c("Tank17"), "Ambient",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(juvtempsHIGH$Tank %in% c("Tank19"), "High",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank20"), "High",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank21"), "Ambient",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank22"), "High",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank23"), "Ambient",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank24"), "High",NA))))))))

juvtempsHURR<-subset(juvtemps, Period==c("Hurricane"), c(Date, TimeDay, Phase, Tank, Temperature, Period))
juvtempsHURR$Treatment<-ifelse(juvtempsHURR$Tank %in% c("Tank17"), "Ambient",
                            ifelse(juvtempsHURR$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(juvtempsHURR$Tank %in% c("Tank19"), "High",
                            ifelse(juvtempsHURR$Tank %in% c("Tank20"), "High",
                            ifelse(juvtempsHURR$Tank %in% c("Tank21"), "Ambient",
                            ifelse(juvtempsHURR$Tank %in% c("Tank22"), "High",
                            ifelse(juvtempsHURR$Tank %in% c("Tank23"), "Ambient",
                            ifelse(juvtempsHURR$Tank %in% c("Tank24"), "High",NA))))))))

juvtempsGROW<-subset(juvtemps, Period==c("Grow-Out"), c(Date, TimeDay, Phase, Tank, Temperature, Period))
juvtempsGROW$Treatment<-ifelse(juvtempsGROW$Tank %in% c("Tank17"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(juvtempsGROW$Tank %in% c("Tank19"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank20"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank21"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank22"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank23"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank24"), "Ambient",NA))))))))

juvtemps<-rbind(juvtempsHIGH, juvtempsGROW)
juvtemps<-rbind(juvtemps, juvtempsHURR)

TimeSeries1j<-ggplot(juvtemps, aes(x = Date, y = Temperature)) + 
  geom_line(aes(color = Treatment), size = 1) +
  scale_color_manual(values=c("blue", "red")) +
  ylab("Temperature (°C)")+
  xlab("Date")+
  ylim(26,32)+
  geom_vline(xintercept = as.POSIXct(as.Date(c("2018-07-31 10:00:00"))), linetype="solid", 
                color = "black", size=1)+
  geom_vline(xintercept = as.POSIXct(as.Date(c("2018-08-16 00:00:00"))), linetype="dashed", 
                color = "black", size=1)+
  geom_vline(xintercept = as.POSIXct(as.Date(c("2018-09-20 16:00:00"))), linetype="dotted", 
                color = "black", size=1)+
  theme_classic()+
  scale_x_datetime(limits = as.POSIXct(as.Date(c("2018-07-30 09:00:00", "2018-09-20 21:00:00")))) +
  theme(legend.position="none")+
  theme(text = element_text(size = 18, color="black"))+
  theme(axis.text = element_text(size=16, color="black"));TimeSeries1j

Plot by time of day by summarizing the mean temp at each time of day for each tank. Display the mean temperatures of each treatment during the Stress test phase (excluding the pause in treatment due to the hurricane system shutdown).

#subset juv temps for the stress period
juvtempsSTRESS<-juvtemps
juvtempsSTRESS$Period<-ifelse(juvtempsSTRESS$Phase %in% c("Stress"), "Stress",
                            ifelse(juvtempsSTRESS$Phase %in% c("Stress2"), "Stress", 
                            ifelse(juvtempsSTRESS$Phase %in% c("Acclimation"), "Acclimation",
                            ifelse(juvtempsSTRESS$Phase %in% c("Hurricane"), "Hurricane",
                            ifelse(juvtempsSTRESS$Phase %in% c("Ambient"), "Grow-Out",NA)))))

juvtempsSTRESS

juvtempsSTRESS<-juvtempsSTRESS[which(juvtempsSTRESS$Phase == "Stress"),]


QuarterlyJ <- ddply(juvtempsSTRESS, c("TimeDay", "Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Temperature, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )

MeansJtemp <- ddply(juvtemps, c("Period", "Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N)
                   )
MeansJtemp

MeansJtemp2 <- ddply(juvtemps, c("Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   max = max(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N)
                   )
MeansJtemp2

Generate final figure.

JuvTempsDay<-ggplot(QuarterlyJ, aes(x = TimeDay, y = mean, color=Treatment, fill=Treatment)) + 
  geom_line(color="black", size = 1) +
  #facet_wrap(~Period) +
  scale_x_continuous(breaks=c(0, 720, 1440), labels=c("0:00", "12:00", "24:00")) +
  scale_color_manual(values = c("blue", "red")) +
  scale_fill_manual(values = c("blue", "red")) +
  theme_classic() +
  geom_ribbon(aes(ymin=QuarterlyJ$lower, ymax=QuarterlyJ$upper), linetype=2, alpha=0.4, color=NA) +
  ylab("Temperature (°C)") +
  xlab("Time of Day") +
  ylim(26,32)+
  theme(text = element_text(size = 18, color="black"))+
  theme(panel.margin = unit(2, "lines"))+
  theme(axis.text = element_text(size=16, color="black")); JuvTempsDay

Combine figures for temperature.

temp_figs<-plot_grid(TimeSeries1j, JuvTempsDay, labels = c("A", "B"), ncol=2, nrow=1, rel_heights= c(1,1), rel_widths = c(1,0.8), label_size = 20, label_y=1, align="h");temp_figs

ggsave(filename="Figures/temp_figs.png", plot=temp_figs, dpi=500, width=14, height=6, units="in")